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A MONTE CARLO PROGRAM FOR CALCULATION OF DOPPLER COEFFICIENTS 


by Robert E. Sullivan 
Lewis Research Center 

SUMMARY 

A Monte Carlo program for the calculation of Doppler coefficients in arbitrary neu- 
tron flux spectra is presented. The program is designed to determine the effective cross 
sections of resonance scatterers and absorbers with a cylindrical or infinite slab geome- 
try. The code may also be used to study the effect of resonance overlap in mixtures of 
isotopes and to investigate the approximations used in analytical calculations of the effec- 
tive resonance integrals. As a test of the accuracy of the program, a flat in lethargy 
spectrum is used to obtain Doppler coefficients for several gold 197 samples, and the 
results are compared with values obtained by using the Nordheim calculation, which is 
expected to be nearly exact for such cases. 

These calculations are compared with published experimental values of Doppler co- 
efficients for gold 197. While a large discrepancy is observed between the Monte Carlo 
and experimental results, the Monte Carlo values agree well with the results of the 
Nordheim calculation. 


INTRODUCTION 

A Monte Carlo program is described that was developed to calculate reaction rates 
and Doppler coefficients for isotopic mixtures of resonance scatterers and absorbers 
with a cylindrical or infinite slab geometry. The basic quantity obtained from the Monte 
Carlo calculations is an ’’effective" cross section which, when multiplied by the inte- 
grated flux existing in the absence of the resonances, yields the correct reaction rate for 
the absorber. The ability to calculate this cross section in any given neutron flux spec- 
trum is of importance in the study of temperature-dependent reactivity effects of non- 
thermal reactors. The code also permits the determination of effective cadmium cutoff 
energies and various activation ratios, as measured by threshold detectors. 

The program evolved from a study of the application of Monte Carlo methods to the 
problem of detecting relatively small changes in reaction rate due to temperature varia- 



tions. From the results obtained, it may be concluded that the method is capable of iso- 
lating these changes accurately. 

The Monte Carlo method is not to be undertaken lightly since it requires much longer 
machine computation times than would an analytical program. Its use is justified only in 
situations where no complete analytical solution is available. Some of the factors affect- 
ing computation times are discussed briefly in the CONCLUSIONS section. 

There are existing Monte Carlo codes which calculate resonance escape probabilities 
and capture rates, but these codes assume an asymptotic solution for the neutron slowing 
down density (refs. 1 and 2). 

The present program was designed to solve the space- and energy-dependent problem 
for simple geometries while making as few approximations as possible. It allows the use 
of an arbitrary input neutron flux spectrum incident on, or distributed through, a sample 
consisting of up to five absorbing isotopes. The individual isotopic contributions to the 
reaction rate are separated so that resonance overlap effects may be evaluated. The 
program operation affords an accurate method for determining differences in reaction 
rate in parallel runjs, as between two temperatures or for various isotropic mixtures, by 
forcing neutrons to follow identical histories to that point where the cases differ. 

The type of Doppler experiment usually performed consists of a temperature- 
dependent activation or reactivity measurement on an absorber material inserted into a 
test region in which the unperturbed flux, or the flux exterior to the region, is known or 
may be calculated. Much of the experimental data available at present consists of mea- 
surements of effective resonance integrals and temperature coefficients in a neutron flux 
which is constant with lethargy in the absence of the absorber. The present Monte Carlo 
code has been used to calculate some of these experiments by generating this neutron 
flux at a boundary removed from the absorber test region. 

There are analytical programs available which will calculate effective resonance 
integrals although they consider the resonances of mixtures of isotopes to be nonover- 
lapping (refs. 3 and 4), or consider the case of overlapping resonances for a 1/E neu- 
tron spectrum only (ref. 5). (All symbols are defined in appendix A. ) 

A test of the precision of the Monte Carlo code is the calculation of Doppler coeffi- 
cients for an isotope such as gold, for which the analytic approximations should be valid. 
This test allows the Monte Carlo results to be compared with calculated results of the 
other programs (refs. 3 and 4) as well as with experiment. 

The analysis section defines the basic quantities required and discusses some of the 
inherent differences in the analytical and experimental definitions of the fluxes used. A 
description of the manner in which individual neutron histories are traced is included 
along with the equations used in selecting neutron paths at points of interaction. The 
equations used in defining capture rates and resonance integrals in terms of the quanti- 
ties computed by the code are also given. 
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The Monte Carlo results obtained for the Doppler coefficients of gold 197 are com- 
pared with values obtained by use of two of the analytical codes (refs. 3 and 4) currently 
in use as well as with experimental results. 

Appendix B describes the input and output of the Monte Carlo program, outlines the 
method used in running problems which have the same geometry but varying cross sec- 
tions, and demonstrates the manner in which the code may be used to solve certain 
classes of problems. Appendix C contains brief descriptions of the input and output for 
the programs used to generate and Doppler broaden nuclear cross sections for the Monte 
Carlo program. The method used in determining an optimum energy lattice for a given 
isotope is referenced, and listings of all codes are given in appendix D. 


ANALYSIS 

Definitions 

In an activation calculation, the quantity desired is the number of neutron captures 
in each absorber isotope over the energy range specified for the problem. The reaction 
rate for a given process is the product of flux and cross section at all energies over the 
volume of interest and is given by 


R = N a f f VpS*’ V)<r(E)dE dV (1) 

V E 

where is the actual flux existing in the sample. The effective cross section for an 
absorber occupying a finite volume is defined as 


'eff 


f y'<p A (E,V)u(E)dE dV 
= V E A 

^^ M < E - V > dEdV 


( 2 ) 


where the flux in the denominator is that which exists in the absence of the absorber. 
The correct reaction rate is then given by 


R " N a°eff ^^'’’M< E - V > dEdV 


(3) 
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The temperature coefficient may be defined in terms of the differences in reaction 
rate occasioned by temperature changes 


a(T } _ 2[R(T 2 ) - R(Ti )] 

(T 2 - T 1 )[R(T 2 ) + R(T X )] 


(4) 


where the mean temperature T m at which the coefficient obtains is taken to be 


T m = T i + ~ ( T 9 " T i ) 


1 ^ 2 VA 2 


(5) 


These definitions may include, as a limiting case, the use of a neutron flux which is 
obtained from an asymptotic solution to the slowing down density. The flux in the ab- 
sence of the absorber is then 


and equation (1) becomes 


<P M (E) = 


q 

?Vs E 



( 6 ) 


The integral in this expression now corresponds to the usual definition of the "resonance 
integral" (ref. 6). This definition requires that the 1/E spectrum not be distorted by 
the absorber, which implies extreme dilution of the absorber. For noninfinitely dilute 
absorbers, the resonance integral must be redefined to allow for arbitrary deviations 
from the assumed 1/E flux spectrum. The effective resonance integral is, from refer- 
ence 7, "the lethargy integral of that cross-section which, when multiplied by the flux 
which would exist in the absence of the resonance, would give the true absorption rate. " 
The effective cross section defined by equation (2), modified to include the neutron 
spectrum of equation (6), then yields the effective resonance integral when integrated 
over lethargy 


X eff - °eff 



f yV A (E,V)a(E)dE dV 
V E 


( 8 ) 
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A given experiment for the determination of the effective resonance integral gener- 
ally consists of placing the absorber, surrounded by a cadmium cover, into a moderator 
region in which the neutron flux, in the absence of the absorber, is approximately uni- 
form in space and lethargy. Experimentally, the flux in the absence of the resonances 
cannot be measured, and the flux must be that which exists when only the absorber poten- 
tial scattering is used or when the absorber is replaced with moderator. Either of these 
fluxes may be used in the present program. The flux incident on the absorber may be 
presumed to be inversely proportional to the neutron energy, but this is not strictly cor- 
rect as the absorber may perturb the flux exterior to its surface (ref. 8). Analytically, 
the limits on the lethargy integral, as defined previously, are taken to be from 0. 4 elec- 
tron volt (0. 64x10"^ J) (the nominal cadmium cutoff energy) to infinity. However, as 
the position of the cadmium cutoff may affect the value of the integral, it is necessary to 
determine this energy accurately for each configuration. The infinite upper limit may be 
replaced by an energy beyond which the absorption rate increases by a negligible amount. 

K equation (8) is used to obtain the effective resonance integral of an absorber at two 
temperatures, the temperature coefficient for that absorber may be calculated in terms 
of the I ^ for suitably small increments in absorber temperature T as 



since 

R J' du 

I e ££ = - — - = Constant x R (10) 

N a^-£'?M< E > v > dEdv 

This temperature coefficient comprises two major effects, the Doppler broadening due 
to a variation in relative velocity of the neutron and target nucleus and that occasioned by 
a change in the geometric size of the sample. The relative importance of the two effects 
may vary greatly and is dependent on the resonance structure and flux spectrum used. 

In the determination of effective resonance integrals for the isotope studied in this re- 
port, where a 1/E neutron flux spectrum is specified, only that component of the tem- 
perature coefficient due to Doppler broadening is important, and the temperature coeffi- 
cient may be replaced by the Doppler coefficient, as defined in equation (9). 
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Monte Carlo Procedure 


The Monte Carlo procedure may be applied so that resonance integrals and, conse- 
quently, Doppler coefficients may be calculated. The geometry and composition of ab- 
sorber selected may be related to specific experiments. It is assumed that there exists 
a test region, which normally contains moderator material, at the outer boundary of 
which the neutron flux is isotropic and uniform in space and lethargy. An absorber sam- 
ple, clad in cadmium is placed into this test region. The test region is large enough so 
that the flux at the outer boundary is not perturbed. Neutrons from the incident distribu- 
tion are then traced through the test region until they are absorbed, moderated, or until 
they escape. Tabulations by energy interval and for each isotope in the absorber are 
made of absorptions and reaction rates for all neutrons entering the absorber. The sum 
of the reaction rates for neutrons of all energies within the absorber volume is the total 
reaction rate. This quantity, divided by the absorber atomic density and the integrated 
flux in the absence of the absorber, multiplied by the lethargy integral, is the effective 
resonance integral, as defined in equation (8). The flux in the absence of the absorber is 
determined by running the program with an identical geometry but with moderator or ab- 
sorber potential scattering replacing the absorber cross sections. 

Although the Monte Carlo code is adaptable to other geometries, the version pre- 
sented in this report is restricted to concentric, finite cylinders within the test region. 
The test region is a surrounding region with planar boundaries, generally containing 
moderator material, and furnishes the surface on which the flux is prescribed, hi order 
to analyze samples whose surface -to-mass ratios approach that of the infinitely dilute 
case, the cylindrical geometry must approach the shape of either a wire or a thin foil. 

In either of these cases, the statistical variance may be greatly increased due to large 
contributions to the reaction rate made by neutrons traveling in the direction parallel to 
the larger cylinder dimension. As most experiments use thin foils, a provision has been 
made in the present program to run the foil and investigate only that component of the 
results perpendicular to the foil surface. 


Neutron Tracking 

A brief outline of the methods and equations used in generating and tracking neutrons 
is presented. The subroutines are described in the order normally encountered in the 
execution of the program. 

INTPOS. - The initial position of the neutron is first specified. As the neutron spa- 
tial density is uniform, the relative frequency with which neutrons appear on a particular 
face is proportional to the ratio of the area of that face to the total area of the outer 


6 



boundary surface. By comparing area ratios with a random number, an initial face is 
chosen. The remaining coordinates are then determined depending on the geometry. For 
a rectangle, for example, if (51 is a random number, 


A + A A x + A v + A z 

-5 Z «ft <_* Z 5 


places the neutron on the Z surface. The x- and y-coordinates are then chosen by ran- 
dom number 


x = X(1 - 2<ft) 
y = Y(1 - 2 (ft) 


Each coordinate is allowed to assume any value over its full range (i. e. , -X < x < X). 
By altering the conditions on the position coordinates, the neutron flux may easily be 
placed throughout the test region rather than restrained to the boundaries. 

E START . - The initial neutron energy is now determined. The number of neutrons 
that has a particular energy is inversely proportional to that energy. The equation used 
is such that each random number returns an energy with the proper relative frequency. 


E 


s 


= E^ exp ((ft In 


E l 


and, in the limits 0 < (ft < 1, the allowed energy range is that specified for the problem. 
Any arbitrary energy spectrum may be used by changing the preceding equation. 

ISOTRO. - Direction cosines for the neutron are assumed to be isotropic in the labo- 
ratory system for the initial flight. In order to avoid tracking neutrons beyond the outer 
boundary planes, the initial direction cosines are chosen to force the neutron into the test 
region. For example, for an initial position on the positive Z-face, 

"z<° 


indicates motion toward the cylinder. 

CMLAB . - All collisions subsequent to the first are isotropic in the center-of-mass 
system. This subroutine uses the standard equations (ref. 9) to find the energy and di- 
rection cosines in the laboratory system based on isotropic scattering in the center-of- 
mass system. 
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SIGM. - The neutron may travel in any of several regions. In the code, all regions 
other than that occupied by the absorbers are considered to contain only one material. 
This subroutine finds the cross sections in any of the outer regions by interpolating be- 
tween input table values. 

SIGA. - In the absorber region, several isotopes may be present. Subroutine SIGA 
finds the cross sections for each isotope in the same manner as SIGM and also computes 
the total microscopic cross sections for the absorber. 

ISQPIK . - In the absorber region, the isotope with which a neutron interacts is de- 
termined by the ratio of the macroscopic cross section of that isotope to the total macro- 
scopic cross section. If 


(ft < 






the interaction is with the isotope I. 

MFP. - The distance D to be traveled by a neutron in the absence of any interac- 
tion with a boundary is then selected by a random choice of the number of mean free 
paths in the medium: 


D = 


-In (ft 


S 


T 


CYLTCK. - The shortest distance to a geometric interaction is determined by this 
subroutine. The equations governing intersection of a directed line with a cylindrical 
surface (ref. 10) are used to determine whether a boundary is crossed before the nuclear 
collision occurs. If so, the direction cosines and energy of the neutron are unchanged 
but a new distance to nuclear interaction is chosen for the material occupying the new 
region. If the nuclear interaction takes precedence, the choice between scattering and 
absorption is made by comparing 



If true, the neutron is scattered; if not, it is absorbed. 

XOYOZO. - When the shortest distance to interaction is determined, the neutron is 
positioned at the correct point by use of the equations 
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x 2 =x l + uD 

y 2 = Yl + w y D 

z 2 ~ Z 1 + “z D 

and tested to determine if it has escaped from the test region. If not, the entire process 
is repeated until the neutron is absorbed, moderated, or until it escapes. 


DOPPLER COEFFICIENTS FOR GOLD 197 

In order to test the precision of the Monte Carlo program in calculating the rela- 
tively small differences in reaction rate caused by the Doppler effect, a flat in lethargy 
neutron flux spectrum has been used to evaluate the Doppler coefficients of several gold 
samples ranging in square root of surface to mass ratio Vs/M from 1 to 50. These 
coefficients were calculated by using equations (8) and (9). 

The use of an absorber consisting of a single isotope in a flux which is inversely 
proportional to the neutron energy allows the Monte Carlo results to be compared with 
values obtained from the Nordheim calculation (ref. 3), which is expected to be nearly 
exact in such cases. In addition, some preliminary experimental temperature coeffi- 
cients (refs. 11 and 12) are available for comparison. The resonance parameters used 
in the analysis were those measured by Desjardins, et al. (ref. 13) with the BWSL 
(Breit-Wigner Single Level) and DBCS (Doppler Broadened Cross Sections) programs 
(described in appendix C) used to generate Doppler broadened cross sections for the 
Monte Carlo code. The experiments referenced were performed on gold foils at temper- 
atures between 295° and 575° K, while the calculations were for foils at 295° and 543° K. 
As the Doppler coefficient is a slowly varying function of temperature, this comparison 
should be valid. The measured reaction rates were used directly, in accord with equa- 
tion (4) 


2[R(T 2 ) - R(Tj)] 

D ~ (T 2 - Tj)[R(T 2 ) + R(Tj)] 

to evaluate the Doppler coefficient and, therefore, no direct value of the effective reso- 
nance integrals for these measurements is given. The experimental error is quoted at 
about 12 percent. 

The results for all methods are given numerically in table I, which presents calcu- 
lated and experimental values of the Doppler coefficient as a function of the Vs/M. As 
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the Doppler coefficient depends on the normalizing value of the effective resonance inte- 
gral, consistent values, as measured by the author of the experimental Doppler coeffi- 
cients, of the effective resonance integral (ref. 14) were used. These values are also 
presented in table I. 

Figure 1 presents curves of the Doppler coefficient as a function of the \/S/M. The 
solid line represents values obtained by the Nordheim calculation with the Monte Carlo 
and experimental values also shown. 

The agreement with experiment, as demonstrated in this figure, is poor. Phenome- 
nologically, it may be expected that the Doppler coefficient curve for gold should have a 
peak, at some relatively high value of the -y/S/M, caused by self -shielding of the large 
first gold resonance, and then rise gradually as the effect of remaining higher energy 
resonances becomes predominant. The experimental measurements suggest only a 
slight variation with surface to mass ratio and agree reasonably well with the Monte 
Carlo results only at lower yS/M values. In an attempt to determine whether the dis- 
agreement was in the Monte Carlo calculation, the quantities AI g ^ and a D were calcu- 
lated by using the analytical code ZUT cited in reference 3. The Nordheim calculation, 
on which this code is based, is expected to give accurate results under the stated as- 
sumptions. Values of the Doppler coefficient calculated by the analytical code TRIX 
(ref. 12) were also compared with the Monte Carlo and ZUT results. As the ZUT and 
TRIX codes use similar analytic approximations, their results for should be consis- 
tent. They do not, however, agree as the TRIX results are consistently lower than the 
ZUT values which are intermediate to, but in better agreement with, the Monte Carlo 
calculations. In addition, although the peak in the Doppler coefficient curve occurs at 
roughly the same surface to mass ratio for all methods, the TRIX peak in AI^ is at a 
higher \/S/M ratio than indicated by the Monte Carlo and ZUT results. 

Nearly all the discrepancy between the Monte Carlo and the Nordheim method (ZUT) 

is attributable to the small differences in Doppler broadened input cross sections used. 

The ZUT program, when run at a surface to mass ratio corresponding to infinite dilution, 

yields a Doppler coefficient of approximately zero. The codes used to prepare the Monte 

-24 2 

Carlo cross sections show an increase of 2 barns (2x10 cm ) in the infinitely dilute 
resonance integral, nearly all of which occurs in the 4. 906-electron-volt (7. 86x10 -J) 

resonance. While the infinitely dilute values would be expected to vary with temperature 
due to weighting with a nonsymmetrical 1/E flux, the magnitude of the change is small 
and depends strongly on the energy intervals used in the numerical integration but infers 
a nonzero infinitely dilute Doppler coefficient. A previously reported change of 0. 7 barn 
(0. 7x10”^ cra^) in 1580 barns (1. 58x10”^ cm^) for I ^ of gold between absolute zero 
and room temperature has been calculated elsewhere (ref. 15). This increase in the in- 
finitely dilute integral is of importance as it is retained for the thicker samples and is 
one of the causes for the higher Doppler values obtained by the Monte Carlo code. In 
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order to verify this, the cross sections used in ZUT for the gold resonance at 4. 906 elec- 
tron volts (7. 86x10" J) were run in the Monte Carlo code at surface to mass ratios for 
which the entire Doppler coefficient is a result of this resonance. Comparison of these 
results, at the -^S/M values indicated in table n, shows that both effective resonance 
integrals and Doppler coefficients are in close agreement. As the effect of this reso- 
nance is less pronounced in the thicker samples, the Monte Carlo and ZUT results, as 
indicated in figure 1 and table I, are in better agreement for the other surface to mass 
ratios investigated even though the cross-section sets used in these two codes differ 
slightly. The results of the Monte Carlo program conform to the shape of the Doppler 
coefficient curve predicted by the analytical code and agree to within about 10 percent 
with the predicted numerical values. 

The differences between the calculated and experimental results are more difficult to 
explain. They are not a result of incorrect cross sections since the most pronounced dis- 
crepancies occur at surface to mass ratios where the Doppler coefficient is dominated by 
the large first resonance for which the parameters are well documented. 


CONCLUSIONS 

The Monte Carlo method outlined in this report is capable of accurately calculating 
small differences in the reaction rate of absorptive samples resulting from a given 
change in temperature. This was demonstrated by using a flat in lethargy neutron flux 
spectrum to calculate Doppler coefficients for a wide range of gold samples. The use of 
this flux spectrum and an absorber which conforms to the approximations required by the 
Nordheim calculation has allowed the Monte Carlo results for the Doppler coefficients to 
be compared with those obtained from an existing analytical code (ZUT) which is known to 
yield valid results. 

The published experimental Doppler measurements exhibit little agreement with 
either the ZUT or Monte Carlo results. The absence of any variation in the experimental 
Doppler coefficient, particularly at surface to mass ratios for which a maximum in this 
function is predicted by both the Monte Carlo and Nordheim calculations, would tend to 
indicate that the experimental measurements should be reevaluated. 

The major disadvantage of the Monte Carlo method is, of course, the large amount 
of machine time required. The running time is a function of the number of resonance 
levels which must be statistically sampled for a given absorber surface to mass ratio 
and of the probable error allowed. While integral quantities, such as reaction rates, 
may be obtained in from 5 to 15 minutes for the range of surface to mass ratios investi- 
gated in this report, the need for extreme accuracy in determining temperature effects 
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requires runs of the order of 15 to 60 minutes, for the determination of Doppler coeffi- 
cients on an IBM 7094 computer. 

The primary advantage of the Monte Carlo method is its flexibility in studying effects 
such as resonance overlap in mixtures of isotopes. The complete solution of the space 
and energy dependent problem in an arbitrary neutron spectrum does not appear to be 
feasible by any other method. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, July 20, 1967, 

120-27-06-18-22. 
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APPENDIX A 


SYMBOLS 


A 

area, cm 2 

X,Y,Z 

initial (surface) coordi- 

Aip 

o 

total surface area, cm 


nates, cm 

D 

distance to any interaction, 

x,y,z 

Cartesian coordinates, cm 


cm 

a 

temperature coefficient, 

E 

energy, eV; J 


1 

rH 

1 

W 

o 

'w' 

E h 

upper energy limit, eV; J 

a D 

Doppler coefficient, (°K) - 

E, 

lower energy limit, eV; J 

h 

average logarithmic decre- 

L 

E s 

*eff 

initial energy, eV; J 
effective resonance inte- 

S s 

ment 

macroscopic scattering 
cross section, cm - * 

AI eff 

gral, b; cm 

change in effective reso- 

S T 

macroscopic total cross 
section for region, cm - * 

N 

nance integral 
total number of neutrons in 

vi 

S T 

macroscopic total cross 
section for isotope, cm - * 

N a 

problem 

absorber atomic density, 

o 

cr 

microscopic cross section, 
cm 2 

N s 

atoms/cm 

moderator atomic density, 

°eff 

effective cross section, 
cm 2 

P.E. 

atoms/cm 2 
probable error 

CT s 

microscopic scattering 
cross section, cm 2 

q 

slowing down density, 

3 


neutron flux in absorber. 

R 

neutrons/(cm )(sec) 
reactions 

neutrons/ (cm ) (sec) 
neutron flux in moderator, 

(R 

random number 

neutrons/ (cm 2 )(sec) 

T 

temperature, °K 


direction cosines for x-, 

TL 

track length, cm 

y-, and z-directions 


T m mean temperature, °K 


u lethargy 

3 

V volume, cm 
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APPENDIX B 


MONTE CARLO CODE INPUT AND OUTPUT 

The input and output for the Monte Carlo program is described. These may vary 
slightly for different versions of the code and the order listed corresponds to the version 
in appendix D. 


Input 


Card 

Quantity 

Format 

Card 

columns 

Remarks 

1 

EH 

E12.5 

1-12 

Upper energy limit 


EL 

E12. 5 

13-24 

Lower energy limit 


XCAP 

E12.5 

25-36 

Half dimensions of outer boundary in x-, y-, and 


YCAP 

ZCAP 

E12.5 

E12.5 

37-48 

49-60 

z -directions 

2 

NHIST 

16 

1-6 

Number of neutron histories 


MIST 

16 

7-12 

Internal mechanism for running more than 32 767 
histories; use any integer less than this 


In each standard case, the number of outer loops run is printed out. In subsequent paral- 
lel runs, MIST is not used directly but must have a value such that NHIST/MIST equals 
the number of outer loops. 


NREG 

16 

13-18 

Number of zones 

NISO 

16 

19-24 

Number of elements or isotopes in absorber 

IGRP 

16 

25-30 

To minimize search time in cross section tables, 

ii. 

only each IGRP energy is tested until neutron 
energy is bracketed 

INTER 

16 

31-36 

Number of neutron histories processed between 
time checks 

NG 

16 

37-42 

Number of energy groups 

PHIBAS 

E12.5 

43-54 

Flux in absence of absorber; may be determined 
in prior run with moderator only present 
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Card 

Quantity 

Format 

Card 

Remarks 




columns 


2 

TEND 

E12.5 

55-66 

Machine time; checked every INTER neutrons and 





final edit performed when actual machine time 
exceeds TEND 

3 

RCYL(I) 

E12.5 

1-12 

Radius and half -height of concentric cylinders. 


H(I) 

E12.5 

13-24 

Innermost first with three sets per card. 

4 

N(I) 

E12.5 

1-12 

Final energy lattice number, abundance weighted 


D(I) 

E12.5 

13-24 

nuclear density, and isotopic mass of each 


A (I) 

E12.5 

25-36 

material; one set per card 

5 

EG(I) 

6E12.5 

1-72 

NG+1 energies bounding NG groups in descending 





order 

6 

E(I) 



Beginning with first isotope in absorber and con- 


s(D 



tinuing through last zone, and energies, scatter 


T(I) 



ing, and total cross sections 


The energies and cross sections are normally read in binary form, in descending order. 
In this case, all energies are first read in, then all scattering cross sections, and, last, 
all total cross sections. Internally, each set of E, S, and T is used with no distinction 
made for different isotopes or materials; that is, if the number of the last E, S, and T 
for the first isotope is i (so that N(l) = i), the first set for the next isotope is i + 1. 
The program differentiates between cross-section sets only by the final energy lattice 
number N(I) which must be correct. 

Decimal cards may be used to read in the cross-section sets if the three BCREAD 
statements in the main program are replaced by: 

READ (5,510) (E(I), S(I), T(I), I = M1,NN) 

510 FORMAT(lP3E12. 5) 


7 

MATCH 

13 

1-3 

Parallel run restart control; yes =1; no = 0 


LRN 

012 

4-15 

First random number constant if new case is run; 
usually 000000000001 
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The pseudorandom numbers generated by the program are controlled by an octal constant, 
LRN; that is, if a given LRN is used to begin generating a series of random numbers, 
this series will be duplicated each time that LRN is used as a lead constant. Most tem- 
perature dependent experiments consist of measurements made on a given external con- 
figuration with, theoretically, only the absorber temperature varying. In such cases, the 
random numbers controlling the neutron path would be temperature independent until the 
neutron entered the absorber and then until the neutron encountered a cross-section set 
changed by the Doppler broadening. The program contains an option for storing on tape 
the lead LRN for each neutron which entered the absorber. The use of these numbers to 
'•restart" neutrons in cases which differ only slightly from the original case allows the 
perturbations to be calculated more accurately. The same method may be used when a 
series of absorbers of the same size but different isotopic compositions are compared. 

The LLRN(LAT) are these octal restart numbers. They are written out on tape for 
all runs made with MATCH = 0 and NPCH =1. If MATCH = 1 and the tape is mounted on 
unit 1, the restart case will be run. 


Card Quantity 

Format 


Remarks 

7 NPCH 

13 

16-18 

Tape storage of random number constants; 




store = 1; no tape = 0 

KSLB 

13 

18-21 

Infinite slab control; yes = 1; no = 0 


Output 

The program output is, in general, sufficiently well described by the printed head- 
ings. Some aspects of certain output quantities are discussed to define them more pre- 
cisely and enable the user to evaluate the results better. 

The fundamental quantity evaluated by the program is the neutron flux. The sum of 
the track lengths of neutrons of all energies in a given volume V is the integrated flux 
(ref. 16) 


J f cp( E, V)dE dV = ^ TL 
V E E, V 


(Bl) 


When the integrated flux is normalized by the total number of neutrons N entering the 
volume, the result is an average track per neutron TL in the volume, and it is this 
quantity which is printed out by the program . 

To use equation (2), the integrated reaction rate in V, in accord with the definition 
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in equation (Bl), is taken to be 


E 


TMEjN^cKE) = f f cp( E, V)N cr(E)dE dV 

a. -'y a 


(B2) 


To find the correct ct^, this quantity must be divided by the integrated flux, in the ab- 
sence of the absorber, produced by the same number of neutrons traced in equation (B2). 
To accomplish this, the quantity in equation (B2) is divided by the average track length 
per neutron in the absence of the absorber TLj^ multiplied by the number of neutrons 
run; that is, 


°eff “ 


£ TL(EME) 
N 


NTL M 


(B3) 


This result, when multiplied by the lethargy integral, is the effective resonance integral. 

As the reaction rates and absorptions are separated by isotope in the program, the 
absorptions, the square root of the surface to mass ratio, and effective resonance inte- 
grals for the individual isotopes are returned by the program. It should be remembered 
that these absorptions, reaction rates, and effective integrals for a particular isotope 
are based on the flux existing in the presence of all the isotopes and, thus, may differ 
slightly from the results obtained from a pure sample of that isotope with the same sur- 
face to mass ratio. 

The values printed out for each energy group are the group number, the upper and 
lower energy bounds, the total number of neutrons with an initial energy in that group, 
the total number with that group energy which enter the absorber region, the reaction 
rate per total neutron and per neutron entering the absorber and the effective resonance 
integral to that energy cutoff. This latter quantity is the I eff between the upper energy 
limit and any given group cutoff energy. 

The reaction rates described in the preceding paragraph may be used to determine 
the effective cadmium cutoff energy. The procedure is to run two cases, identical but 
for the presence of cadmium in one. Then, the cutoff energy is that energy in the non- 
cadmium case at which the reaction rate per entrant neutron is equal to the total reaction 
rate per entrant neutron in the cadmium covered case. 

Finally, the sample mean (the effective resonance integral), variance, and probable 
error are printed out. These are based on the individual neutron results and the stand- 
ard equations (ref. 17): 
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Although standard methods may be used in the error analysis for directly sampled 
quantities such as the effective resonance integral, the use of the restart procedure in 
determining difference quantities, such as the Doppler coefficient, complicates the error 
analysis. Although alternate methods are possible, the method used in this report has 
been to set an upper limit to the probable error for restart cases by defining 
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APPENDIX C 


AUXILIARY PROGRAMS 

The auxiliary programs used in selecting an energy - cross-section lattice, gener- 
ating cross sections, and Doppler broadening the cross sections when temperature- 
dependent results are required are described. While any cross-section set with suffi- 
cient detail may be used in the Monte Carlo code, these programs make few approxima- 
tions and should increase the accuracy of the calculation. There are three auxiliary pro- 
grams. 

EPIGRAM chooses an energy lattice spaced so that no more than a specified error is 
involved in interpolating the cross section between the lattice points chosen. The error 
specified may be an absolute difference in barns (cm ) or a percent error in the absorp- 
tion or total cross section. 

BWSL generates cross sections by using the single level Breit-Wigner equation. To 
avoid the incorrect treatment (ref. 18) of interference between resonance and potential 
scattering, the equation is coded in complex number form. 

DBCS Doppler broadens the cross sections obtained from BWSL. As the EPIGRAM 
code returns an energy lattice in which little error is incurred by interpolating between 
cross-section points, the DBCS code uses interpolated cross sections for numerical inte- 
gration rather than require an evaluation of the single level Breit-Wigner equation at 
each point. 

A brief resume of the equations, input and output, and comments on each of these 
programs, with the exception of EPIGRAM which is referenced elsewhere, follows. 


BWSL 

The cross sections used in the main program are generated by using the single level 
Breit-Wigner equations. The complex number forms of the equations were coded to 
avoid the occurrence of negative scattering cross sections which result from an incorrect 
treatment of the resonance and potential scattering interference. The equations for the 
scattering and absorption cross section are 
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Equation (C2) is also used for fission cross sections when r, is replaced by . The 
symbols used denote 


A target mass/neutron mass 
B isotopic abundance 
c capture 

E incident neutron energy 
f fission 

g statistical weight factor 

J summation over all resonances in an isotope having same g k 
k summation over all statistical weight factors 
o value at resonance 
R nuclear radius of target isotope 
s scatter 

6 phase shift, taken to be -R/ft 
r nuclear level width 
ft rationalized wavelength 
a microscopic cross section 

The scattering cross-section equation (eq. (Cl)), due to summing the resonance part 
over all statistical weight factors with an attendant error in the potential scattering con- 
tribution, cannot be used directly. In the program, the potential scattering multiplied by 
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the statistical weight factor is subtracted from the expression given leaving the resonance 
and interference parts. When the summations are completed, the total potential scatter- 
ing for that isotope is then added. 


NOGS 

Input 

number of different g factors in each isotope 

AA 

atomic mass of each isotope 

ABB 

abundance of each isotope 

RR 

nuclear radius of each isotope 

LEND 

number of isotopes 

LEND 

number of energy points 

E 

energies 

G1 

g value of set of resonances 

JEND 

number of resonances with this g value 

ERI 

resonance energies 

GAMG1 

radiation widths 

GAMN1 

neutron widths 

GAM FI 

fission widths 


DBCS 


The Doppler broadening of the cross sections generated by BWSL is accomplished by 
this program. Rather than generate cross sections with the BWSL code at each point re- 
quired in the numerical integration (eq. (C3)), the zero degree cross sections are inter- 
polated, and a Simpson's rule integration with up to 400 points is performed. The equa- 
tion used assumes a Maxwellian distribution for the target nuclei, and the integration is 
over all energies within which the exponential term e -x lies within the limits 
-88 < x < 88 
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where k is the Boltzmann constant, T is the temperature in °K, E n is the neutron 
energy at which the temperature dependent cross section is to be found, K is a constant, 
and the remaining symbols are the same as those defined for the BWSL code. 

In order to obtain an accurate value for the infinitely dilute integral, the cross sec- 
tions are assumed to vary linearly between energy points and the 1/E flux weighted 
resonance absorption integral 



is printed out at each energy interval. 


(C4) 


Input 


A 

T 

TEND 

NLEV 

ITYPES 

NPUNCH 
NEDU 
NED L 


atomic weight 
temperature, °K 

machine time to termination of run and final edit 
number of energy levels 

flag for number of types of cross sections to be broadened; -1, 0, 1 for 
1, 2, or 3 types 

if set equal to 1, the broadened cross sections are punched on cards 
upper energy at which broadening is performed 
lower energy at which the broadening is performed 


As the integration relies on existing cross sections over an input energy range from 
NEDU to NEDL, the first and last few points may be in error. If this occurs, an error 
return states the fact and gives the energy cutoff actually used. 


INTER number of energies processed between time checks 

E energies in descending order 

SIGC capture cross sections 

SIGS scattering cross sections 

SIGF fission cross sections 


The output consists of the energies with the corresponding broadened cross sections 
and resonance absorption integrals. 
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APPENDIX D 


PROGRAM LISTINGS 


Monte Carlo Program ERIDE 


S I B J OB 

& I BFTC ERI LIST 

DIMENSION E I 5000 ),S( 5000) ,T( 5000) * S I ( 8 ) , T I ( 8 ) , RCVL ( 3 ) , H( 3 ) , D ( 8 ) , 

1 A < e ) ,N(8> ,KKK(8 ) ,ARATE(5),RH0(5),SMR(5),RRIN(5),ERI(5) ,NA6S!8) , 

2 EG I 100 ) , EDR ( 100 ) , SEDR ( 100 ) , ER I P ( 100 ) , NEST ( 100) ,NHST ( 100 ) , DELU ( 100 ) 
3,GDIST( 100), LLRN! 7500) 

COMMON E,S,T 
COMMON/ RDEP/LRN 
REWIND 1 

READ ( 5,501 ) EH,EL,XCAP,YCAP,ZCAP 

501 FORMAT! 1P5E12.5 ) 

WRITE(6,500)XCAP,YCAP,ZCAP 

500 FORMAT! 1 HI //// / 10X , 69HTHE OUTER MODERATOR DIMENSIONS, AT WHICH THE 

1 NEUTRON FLUX IS 1/E, ARE/// 14X , 10HX BOUNDARY , 1 OX , 10HY BOUNDARY, 10 

2 X , 1 OHZ BOUNDARY/ //5X »1P3E20«5) 

WRITE !6,508)EH,EL 

508 F0RMAT!1HL,10X,38HTHE ENERGY RANGE FOR THIS CASE IS FROM , 1 PE 1 2 . 5 , 2 
1X,2HT0,1PE12.5,2X,15HELECTR0N VOLTS. ) 

READ! 5,505 )NHIST ,MIST ,NREG,NISO, IGRP, I NT ER , NG , PH I B AS , T END 

505 FORMAT ( 7I6,1P2E12.5 ) 

WRITE ! 6,506 ) NH I ST , NR EG , N I SO , I GRP , I NT ER , PH I B AS , T END 

506 FORMAT! 1 HL , 1 OX , 48 HTHE INPUT DATA USED FOR THIS CASE IS AS FOLLOWS. 

1///116H HISTORIES NUMBER OF NUMBER OF ENERGI 

2ES IN TIME FLUX IN ABSENCE EST I M AT ED/ 1 1 8H 

3RE0UESTED REGIONS ABSORBERS MAJOR GROUP 

4 INTERVAL OF ABSORBER RUNNING T I ME/ / 1 1 1 , 2 1 1 6 , I 1 7 , I 1 7 

5,1PE24.5,1PE18.5 ) 

NUTS-NREG— 1 

NTOT=NISO+NUTS 

NEG-NG+1 

READ (5,502 ) (RCYL ( I ) , H ( I ) , 1=1, NUTS) 

502 FORMAT! 1P6E12. 5 ) 

WRITE (6,503) (I,RCYL(I),H(I) , 1=1, NUTS) 

503 FORMAT! 1 HL , l OX , 49HTHE DIMENSIONS OF INNER CYLINDRICAL REGION NUMBE 
1R, I3,2X,3HARE,///20X,6HRADIUS, 1PE13. 5,2X, 11HHALF HE I GHT , 1 P E 1 3 . 5 ) 

READ ! 5 , 504 ) (N(I),D(I),A(I), I = l , NTOT ) 

504 FORMAT! 112, 1P2E12. 5) 

WRITE(6,507J MI SO,NTO T , ( I , N ( I ) , D ( I ) , A ( I ) , I = 1 , NTOT ) 

507 FORMAT ( 1HL, 10X ,9HTHERE ARE , I 3 , 2X , 46HI SOT OPES IN THE ABSORBER REG 1 0 

IN AND A TOTAL OF , I 3 , 2X , 10HM AT ER I A L S . / / / 62 H MATERIAL FINAL 

2 ENERGY NUCLEAR AT0MIC/62H NUMBER LATTIC 

3E NUMBER DENSITY W E I GHT/ / ( I 9 , I 20 , 1 P E 1 9 . 5 , 1 P E 1 7 . 5 ) ) 

READ! 5,509) ( EG ( J ) , J = 1 , NEG ) 

509 FORMAT! 1P6E12.5) 

M 1 = 1 

DO 5000 J = 1 , NTOT 
NN = N ! J ) 

CALL BCREAD ( E ( Ml ) , E ( NN ) ) 

CALL BCREAD (S(M1 ) , S ( NN ) ) 

CALL BCREAD (T (Ml ) , T ( NN ) ) 

M 1 =NN+ 1 
5000 CONTINUE 

READ! 5,512 ) MATCH , LRN ,NPCH , KSLB 
512 FORMAT! 13,012,213) 

5050 CALL T I ME 1 ( TST ) 

TSTART=TST /3600.0 
CN=ALOG( EH/EL ) 

NM0D=0 
NES P = 0 
NTAB=0 
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KENT=0 
NENT-0 
NEECL=0 
AD IST = 0 *0 
NA B 1 =0 
NAB2=0 
NA B 3 = 0 
NA B4 = 0 
NA B 5 = 0 
XX=0.0 
XS=0.0 

MISTY = MIST 
DO 5001 J=l f NG 
NEST (J) = 0 
NHST(J) = 0 
GD I ST { J )=0.0 
EDR ( J ) =0 .0 
5001 SEDR<J>=0.0 

DO 5003 L = 1 ,N I SO 
5003 A RAT E ( L ) =0 • 0 
BD I S T = 0 . 0 
CD I ST = 0 .0 
DD I ST=0 • 0 
RDIST=0.0 

NOUT =(NHIST/MIST )+l 
NREM=NHIST-MIST*(N0UT-1 ) 
JINTER=INTER+MIST 
DO 1003 10=1, NOUT 
jinter=jinter-mist 
L AT=0 

IFIN0UT.60. IO)GO TO 1007 
IF(MATCH.EQ.O) GO TO 1005 
READ { 1 ) MIST, (LLRN(LAT), LAT=1,MIST) 
GO TO 1005 
1007 MIST = NREM 

I F ( NREM .1 E *0 ) GO TO 1003 
I F ( MATCH • EQ* 0 ) GO TO 1005 
RE AD 1 1 ) MIST,(LLRN(LAT),LAT = 1,MIST) 
1005 DO 1000 1 = 1, MIST 

I F ( I *LT * J INTER )G0 TO 1111 
J INTFR = J INTER+INTER 
CALL TIMEl(TST) 

TPRES=TST/3600.0 

TIME=TPRES-TSTART 

ifitime.lt *tend)go to mi 

NHIST = 110-1 )*MISTY+I-1 
IF(NPCH.EQ.O) GO TO 1101 
KENT =N ENT -KENT 

WRITE (1 ) KENT, ( L IRN < L AT ) , L AT= 1 , KENT ) 
GO TO 1101 
1111 CONTINUE 
KKKU ) = 1 
NK=NT0T-1 
DO 2101 J = 1 , N K 
2101 KKK ( J+l ) =N ( J ) +1 
I NSA P = 0 
K = 0 

M=NREG 
I NRU P= 1 
XV=0.0 



520 I F ( K— 1 >525,535,535 
525 IF(MATCH.E0.0)G0 TO 530 
LAT=LAT+1 
LRN=LLRN ( L AT ) 

530 CALL IN IPOS ( XC AP , YC AP , ZC AP , XO , YO, ZO , IN, IGN,KLRN> 
ELA=ESTART (EL,EH,CN) 

DO 536 IHST=1,NG 

IF(ELA.LT.EGUHST).AND.ELA.GE.EG(IHST + I))G0 TO 537 

536 CONTINUE 

537 NHST ( IHST )=NHST ( IHST ) + 1 

535 CALL ISOTRO ( WX » WYt WZ *K» INt IGN) 

IF(K-1 >550,540,540 

540 CALL CMLAB(WX,WY,WZ,ELA,A( II >,CX,CY,CZ> 

IF(ELA-EL >545,545,551 
545 INRUP=2 

GO TO ( 551 ,950,551 ,551 >, INRUP 

550 K= 1 

551 CONTINUE 
CX = WX 

C Y = W Y 
CZ = WZ 

552 CONTINUE 

I F ( M*GT • 1 ) GO TO 710 
N S = 1 

S I GT = 0 «0 
S I GA = 0 • 0 

610 DO 650 L = 1 , N I SO 
NF=N ( L > 

IF(KKK(L>.GT.NS)GO TO 616 
DO 612 K = NS , NF , I GR P 
KK = K 

IF(E(K)-ELA>614,614,612 
612 CONTINUE 

KK = NS+ I GRP 
614 K = KK- I GR P 
GO TO 618 
616 K=KKK ( L ) 

618 DO 620 K = K » NF 
KK= K 

IF ( E (K ) -EL A >630, 6 20, 620 
620 CONTINUE 
630 K=KK-1 

AM= (S(K)-SIK+1)I/(EIK)-E(K+1)) 

8S = S ( K ) —AM*E ( K > 

SI(L)=AM*ELA+BS 

TM=(T(K)-T(K+1 ) )/(E (K)-E(K+1 ) ) 

BT = T ( K ) -TM#E < K ) 

T I ( L )=TM*ELA+BT 
SIGT=SIGT+D(L ) *T I (L> 

SIGA = SIGA + D<L ) *( T I (L ) — S I ( L > ) 

N S = N F + 1 
KKK ( L ) =K 
650 CONTINUE 

CALL ISOPIKtNISO, II ,D,TI , S I GT ) 

GO TO 553 
710 L=NIS0+M-1 
I I =L 

NS=N( L-l >+l 
NF = N ( L ) 

IF(KKK(L).GT.NS)GO TO 716 


DO 712 K-NS *NF , IGRP 
KK=K 

IF(E(K)-ELA)714t714,712 
712 CONTINUE 
KK=NS+ IGRP 
714 K=KK- IGRP 
GO TO 718 
716 K = KKK ( L ) 

718 DO 720 K^K * NF 
KK=K 

IF( E(K J-ELA ) 730*720,720 
720 CONTINUE 
730 K=KK- 1 

L -N I SO + M- 1 

XM=(S(K)-S(K + 1 ) )/<E(K)-E(K + l ) ) 

B X = S ( K ) —X M * E ( K ) 

SI (L)=XM*ELA+BX 

XMT=(T(K)-T ('< + !) >/(E(K)-E(K+l) ) 

BXT=T(K)-XMT*E<K) 

T I (L)=XMT*ELA+BXT 
KKK(L)=K 
553 CONTINUE 

CALL MFP (T I ( I I ) »DIST*D( I I ) ) 

MM = M 
NFLG = 1 
I SE T = 0 

555 IF ( MM-NREG ) 1510* 1521 *2000 

1510 CALL CYLTCK(-ltIFLG,RCYL(MM ) * D I ST * XO * YO * ZO * WX * WY , WZ ) 

I F ( IFLG.E0.~1 ) HZ=H(MM) 

I SET= I FLG 

NFLG= I ABS ( I FLG ) +1 

M=M+ I AB S ( I FLG ) 

1520 IF(MM-1 >2000*1516*1521 

1521 CALL CYLTCK ( +1 , I FLG * RCYL ( MM-1 ) ,H( MM-1 ) * D I ST * XO * YO f ZO * WX * WY , WZ ) 

I F ( IFLG.E0.~1 ) HZ-H { MM- 1 ) 

ISET=ISET+2*IFLG 
M=M-NFLG*IA6S ( I FLG ) 

1516 I F C (MM-1 ) .NE.O. OR. INSAP. NE.OJGO TO 5556 
I NS A P= 1 

DO 5557 IEST=1,NG 

IF( ELA.LT.EG( IEST ) . AND.ELA.GE.EGt I EST+1 ) ) GO TO 555V 
5557 CONTINUE 

5559 NEST( IEST)=NEST ( IESTI+1 
IF(MATCH.NE.O)GO TO 5556 
L AT=L AT+1 
LLRN ( L AT ) =KLRN 

5 556 CALL XO YOZO ( D I ST , WX * WY * WZ * XCAP * YC A P * ZC AP , XO * YO * ZO * HZ * I SET* INRUP) 
GO TO (5550*5552*5553*5554) * M M 

5550 CONTINUE 
IF(KSLB.EO.O)GO TO 5549 
DIST = DIST ;i *ABS(WZ) 

5549 ADI$T=ADI$T+DIST 
XIND=DIST*S IGA 
XX=XX+XIND 
XV=XV+XIND 
DO 5551 L = 1 *N I SO 

ARATE ( L ) = ARAT E ( L ) +D I ST*D ( L ) * ( T I ( L ) -S I ( L ) ) 

5551 CONTINUE 

DO 5560 I EG = 1 * NG 

IF(ELA.LT.EG( IEG) •AND.ELA.GE.EG( I EG+1 ) )G0 TO 5565 
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5560 CONTINUE 

5565 El)RUEG)=EDR ( I EG ) + D 1 ST*S I G A 
GDISTt IEG ) -GD I ST ( IEGJ+DIST 
GU TO 5558 

5552 BDIST=BDISr+DIST 
GO TO 5558 

5553 CDIST=CDIST+DISr 
GO TO 5558 

5554 DDIST=DDIST+DIST 
5558 RDIST=RDIST+DIST 

GO TO (560,560,960,560) , INRUP 
560 IF(MM-M)552,572,552 
572 CONTINUE 

CALL INTACT (SI ( I I ) »TI ( 1 1 ) t INRUP) 

GO TO ( 535,535,535,970) t INRUP 
950 NMOD=NMOD+ 1 

NENT=NENT+INSAP 
GO TO 999 
960 NESP=NESP+1 

I F < INSAP • EQ. 0 )G0 TO 965 
NEECL=NEECL+I 
965 NENT=NENT+ INSAP 
GO TO 999 

970 IFIM-1 ) 980, 974, 980 

974 GO TO {975,976,977,978,979), II 

975 NAB1=NAB1+1 
GO TO 980 

976 NAB2=NAB2+ 1 
GO TO 980 

977 NAB3=NAB3+1 
GO TO 980 

978 NAB4=N AB4+ 1 
GO TO 980 

979 NAB5=NAB5+ 1 

980 NTAB=NTAB+ 1 
NENT=NENT+INSAP 

999 XS=XS+XV*XV 

1000 CONTINUE 

I F ( NPCH. EQ.O ) GO TO 1003 
KENT =NENT -KENT 

WR I TE I 1 ) KENT , (LLRN(LAT) ,LAT-1,KENT) 

K ENT =NENT 
1003 CONTINUE 
1101 WR1TE16,1124) 

1124 FORMAT { 1 HI ) 

IF(KSLB.NE.O)GU TO 1130 
WR I TE ( 6 , 1 1 2 6 ) 

1126 FORMAT { 1 H L , 1 0 X , 46 H TH I S CASE IS FOR A THREE DIMENSIONAL CYLINDER.) 

GO TO 1001 
1130 WRITE(6,1128) 

1128 FORMAT ( 1HL , 10X, 3 4HT HIS CASE IS FOR AN INFINITE SLAB.) 

1001 WRITE (6,1 100)NHIST ,NENT ,NEECL,NAB1,NAB2,NAB3,NAB4,NAB5,NTAB ,NMOD,N 
1ESP 

1100 FORMAT ( 1HL , 1 OX , 37HRESULT S FOR THIS CASE ARE AS FOLLOWS .//// 1 1 9H HI 
IS TORIES ENTERING ESCAPING NEUTRONS ABSORBtD BY ISOTO 

2PE NUMBER TOTAL MODERATED ESC APED/98H RUN 

3 CYLINDER CYLINDER ONE TWO THREE FOUR 

4 FIVE ABSORBED// 18, 21 II, I 10, 2111, 110,112,113,110,112) 

H=ABS ( H) 

PH I T = CN 
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VSAM=6.2832*RCYL*RCYL*H 

SURC=6.2832*RCYL*(RCYL+2.0*H) 

IF(KSLB.EQ.0)G0 TO 1120 

VSAM=H 

$URC=1 .0 

1120 FN=FLQAT ( NENT ) 

FH= FLOAT (NFIIST ) 

DT = 0*0 

ARTOT = 0 *0 

DO 1150 L=1,NIS0 

RHO( L )=D(L ) -A (L J/0.6025 

DT = DT+D ( L ) 

SMR( L )=SORT (SURC/ (RHO(L )-VSAM ) ) 

RRIN(L>=ARATE (D/FLOAT (NENT) 

ERI(L) = (RRIN(L ) *PH I T ) / < D ( L ) *PH I BAS ) 

RH0T=RH0T+RH0(L) 

ARTOT=ARTOT+ARATE ( L ) 

1150 CONTINUE 

WRITE (6, 12 50) ( L * RHO ( L ) * SMR ( L ) , RR I N < L ) , ER I U> *1=1* NISO) 

1250 FORMAT < 1HL ,99H ABSORBER DENSITY SQUARE R 

100T REACTION RATE PER EFFECTIVE RESONANCE/1 2H NUMBER, 34 

2X,48H0F S/M INCIDENT NEUTRON I NT EGRAL // I 1 0 , 6X , 1P4E 

320.5) 

SUME = 0 • 0 

SDELU*=0 . 0 

FNEST=0 .0 

FNHST=0 • 0 

DO 1160 I G= 1 * NG 

FNEST = FNEST + FLOAT (NEST { IG) ) 

FNHST = FNHST + FLOAT (NHST ( IG) ) 

DELU < I G ) =ALOG ( EG ( IG ) /EG ( IG+1 ) ) 

SDELU=SDELU+DELU( IG) 

SUME=SUME+EDR< IG) 

IF (GDI ST ( IG).EQ.0.0) GO TO 1158 
GDIST ( I G ) -E DR ( IGJ/GDIST ( IG) 

1158 SEDR( IG)=SUME/FNHST 
EDR ( IG)=SUME/FNEST 

1160 ERIP( IG)=SUME*SDELU/(DT*PHIBAS*FNEST) 

WR ITE ( 6, 1 105 ) 

1105 FORMAT! 1HL ,10X,93HTHE REACTION RATE IN THE ABSORBER, PER TOTAL ENT 
1RANT NEUTRON, BY ENERGY GROUP, IS AS FOLLOWS.) 

WRITE (6, 1 106 > ( IGtEGUG) , EG ( IG+1 ) ,NHST ( IG) ,NESI ( IG) ,SEDR( IG) , EDR ( I 
2G ) , ER I P ( IG ) , IG=1 ,NG ) 

1106 FORMAT! 1HL,112H GROUP UPPER LOWER NST ART 

2 NENTER NORMALIZED RATE TO LOWER LIMIT TOTAL IEFF T0/111H 

3 NUMBER ENERGY ENERGY IN GROUP P 

4ER TOTAL N PER ENTRANT N LOWER L I M 1 1 / / ( I 6 , 1 P2 E 1 8 . 5 , 2 1 8 , 1 P 

53E18.5) ) 

AD I ST = ADIST/FLOAT (NENT) 

SMRT=SQRT (SURC/ (RHOT*VSAM> ) 

RT I N = ARTOT /FLOAT (NENT ) 

ERIT=(RTIN*PHIT )/ (DT*PHIBAS) 

WRITE (6,1107) AD I ST 

1107 FORMAT ( 1 HL , 65HTHE INTEGRATED FLUX PER ENTRANT NEUTRON IN THE ABSO 
SRBER REGION IS,1PE15.5) 

WRITE (6,1200) PHIT »SURC?VSAM» RHOT , SMRT , RT I N, ER I T 
1200 FORMAT(lHL,10X,20HT0TALS FOR THE CASE.///115H LETHARGY ABS 

10RBER ABSORBER DENSITY SQUARE ROOT REACTION RA 

2TE PER EFFECTIVE RE SON ANCE /40H INTERVAL SURFACE 

3V0LUME ,24X,45H0F S/M INCIDENT NEUTRON I NTEGRAL/ / 1 P E 13 
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4.5,1P2E15«5,1P2E15«5»1PE17.5,1PE22.5) 

FN = FLOAT{'VENT ) 

D E N = P H I T / ( DT*PHIBAS ) 

XX=XX*UEN/FN 

XS= ( XS*DEN*DEN/FN)-(XX*XX ) 

DELX=0.6745*SQRT (XS/IFN-l.O)) 

WKITE(6»1213) XX,XS,DELX 

1213 FORMAT ( 1HL ,42H SAMPLE VARIANCE PROBABLE/lOH 

1MEAN,26X, 5HERROK// 1PE14.5, 1PE15.5,1PE16.5) 

I F ( M A T C H • N E . 0 ) GO TO 1298 
IF(NPCH.EL).O) GO TO 1299 
WRITE ( 6 ,12 7 5 )NENT 

1275 FORMAT ( 1 HL * 10X , 46H0CT AL CAROS HAVE BEEN PUNCHED WHICH WILL FORCF,I 
16,2X,38HNEUTR0NS TO ENTER THE ABSORBER REGION.) 

WRITE(6tl277> 10 

1277 FORMAT ( 1HL ,10X,38HTHE NUMBER OF OUTER LOOPS COMPLETED I$,I4,1H.) 

GO TO 1299 

1298 WRITE (6,1274)NENT 

1274 FORMAT ( 1HL ,10X,23HTHIS CASE WAS RUN US I NG* I 6, 2 X t 77H0CT AL RESTART N 
1UMBERS WHICH FORCED EACH NEUTRON TO ENTER THE ABSORBER REGION.) 

1299 CONTINUE 

WRITE (6r2U20) LRN 

2020 FORMAT( 1 HL t 1 OX , 57HTHE LAST RANDOM NUMBER COEFFICIENT USED IN THIS 
ISERIES IS ,3X ,012 ) 

WRITE (6,9999 ) < GD I ST ( I G ) , I G= 1 , NG ) 

9999 FORMAT ( 1HL ,1P6E12.5 ) 

REWIND 1 
CALL EXIT 

2000 WR I TE ( 6 , 200 1 ) 

2001 FORMA T ( 20H TROUBLE IN INDICES.) 

STOP 

END 

MBFTC GEN LIST 

SUBROUTINE R ANDM ( K , F ) 

common/rdep/lrn 

LRN=LRN*30517578125 

Z R = 2 .91038304567E-11 

R=FLOAT(LRN)*ZR 

RETURN 

END 

$IBFTC PROSIT LIST 

SUBROUTINE I N I POS ( XC AP , YC AP , ZC AP , XO , YO , ZO , I N , I GN T KLRN ) 
COMMON/RDEP/LRN 
I GN = 1 

AX= YCAP^ZC AP 
AY=XCAP*ZCAP 
A Z = XCA P# YC AP 
AT=AX+AY+AZ 
RX = A X/ AT 
RY=AY/ AT 
KLRN=LRN 

CALL RANDM ( R 1 , 1 ) 

X0 = XCA P* ( 1 .0-2 .0*R1 ) 

CALL RANDM(R2,1> 

Y0=YCAP*(1.0“2.0*R2) 

CALL RANDM ( R3 , 1 ) 

Z0=ZCAP*(1.0-2.O*R3) 

300 CALL RANDM ( R4 , 1 ) 

I F ( R4— RX— RY >310,320,320 
310 IF (R4-RX) 340,330,330 



320 Z0=ZCAP*(Z0/A6S(Z0) > 

IN = 1 

IF( Z0.GT.0.0)G0 TO 350 

IGN=“1 

GO TO 350 

330 YO=YCAP*( YO/ABS(YO) ) 

I N=0 

IF( YO.GT.O.O )G0 TO 350 

IGM=-1 

GO TO 350 

340 X0=XCAP* { X0/A6S I XO ) ) 

I N = — 1 

IF< XO.GT.O.O )G0 TO 350 
I GN = -1 

350 CONTINUE 
RETURN 
END 

$ I 6 FTC EOPEN LIST 

FUNCTION ESTART ( EL ♦ E 1 ,CN ) 

COMMON/ RDEP/LRN 
CALL RANDM ( R30 , 1 ) 

ESTART=EL-EXP (R30-CN) 

RETURN 

END 

SI8FTC COR I EN LIST 

SUBROUTINE I SOTRO ( WX , WY , WZ , K , IN, IGN) 
COMMON/ RDEP/LRN 
100 CALL RAN0M(R7,1) 

WZ = R7 + R7-1 • 0 

SNT = SORT ( 1.0-WZ*WZ ) 

120 CALL RAN DM ( R8 , 1 I 
CALL RAN DM ( R9 , 1 ) 

D=R9+R9~1 .0 

E=R8*R8 

F=E+D*D 

IFIF-1. 0)130, 130,120 
130 WX=SNT*(E+E-F )/F 

WY=SNT*(2 .0*R8*D ) /F 
T = W X * WX + W Y * W Y + W Z * W Z 
T = S 0 R T ( T) 

W X = W X / T 
W Y = W Y / T 
W Z = W Z / T 

IF(K-1 >135,190,190 
135 I F < IN)140,150,160 
140 IF( IGN.GT.OJGO TO 142 
IFIWX)100,190»190 
142 IF(WX) 190,190,100 
150 I F ( I GN • GT • 0 ) GO TO 152 
IF (WY) 100, 190, 190 
152 IF(WY)190,1 90 ,100 
160 I F ( I GN • GT • 0 ) GO TO 162 
IF(WZ)100, 190,190 
162 I F ( WZ I 1 90 , 1 90 ,100 
190 CONTINUE 
RETURN 
END 

$ I BFTC LUTATE LIST 

SUBROUTINE CML A B ( WX , WY , WZ , EL A , A , CX , C Y , CZ ) 
COMMON/ RDEP/LRN 



FL B = EL A 
Z = A/1 .00898 

230 CMU=CX*WX+CY*WY+CZ*WZ 

CKON=SQRT ( 1.0+Z*Z+2.0*Z*CMU) 

W X = (CX+Z^WX )/CK0N 
WY={ CY+Z*WY ) / CKON 
WZ=<CZ+Z*WZI/CK0N 

ELA=(ELB*CKDN*CKON)/( (Z+1.0)*(Z+1.0>> 

CONTINUE 

RETURN 

END 

$ I BFTC FLIGHT LIST 

SUBROUTINE MFP < S IGT , D I ST , D ) 

CGMMON/RDEP/LRN 
CALL RANDM ( RIO * 1 ) 

FNMFP“-AL0G(R10) 

SIGMAC=D*S IGT 
2020 DIST=FNMFP/SI GMAC 
RETURN 
END 

$ I B FTC GOK LIST 

SUBROUTINE CYLTCK ( N , I F LG » RCYL , H , SM , XO f YO * ZO , WX , WY , WZ ) 
COMMON/ RDEP/LKN 
SQRTF ( X )=SQRT (X) 

QUADRA ( 0 ) =A*Q**2+2 .0*B*0+C 
I FL G = 0 
H=ABS( H) 

EPS 1 = 1 .OE-4 
IF(N)2,97,3 

2 A=1.0-WZ**2 
B=XO#WX+YO*WY 
C=X0**2+Y0**2-(RCYL*RCYL ) 

IF (WZ >4,5,6 

4 H = -H 

6 T=(-ZO+H)/WZ 
IF(T>97,8,9 

8 I FLG=-1 
GO TO 40 

9 I F ( QUADRA ( T ) )10,10,5 

5 IF ( QUADRA (SM) ) 9 9, 9 9, 11 
11 e=B**2-A*C 

IF(E)12,13,14 

14 SM= (-B + SQKTF ( E ) )/A 
I F L G = 1 

GO TO 99 

10 IF( T-SM ) 15 ,99,99 

15 SM = T 
IFLG=-1 
GO TO 99 

3 IF( ZO+H) 16,16,17 
17 IF(Z0-H)18,19,19 

19 IF(WZ)220,99,99 

16 I F ( WZ ) 99 ,99,20 
220 H=— H 

20 C=X0**2+Y0**2-(RCYL*RCYL ) 

B=XO*WX+YO*WY 

IF(C)21,22,22 

22 IF(B)23,99,99 

21 A=1.0-WZ**2 
U=-( ZO+H) /WZ 



I F ( QUADRA ( U ) >24,99,99 

23 A=1.0-WZ**2 
e=b ** 2 -A*C 

I F { E ) 99 , 99 * 2 5 

25 U=- { ZO+H ) / WZ 

I F ( QUADRA { U ) >24,24*26 

26 IF< A*U+B)27,99,99 

24 IF ( U-SM >28 ,99,99 

28 SM=U 
H~— H 

iflg=-i 

GO TO 99 

18 C=X0**2+Y0**2- ( RCVL-KCYL > 

B=XO*WX+YO*WY 

IF(C)29,30,30 

30 I F ( B ) 31 ,99,99 

31 A=1.0-WZ**2 
E=B--2-A^C 

I F ( E ) 99 ,99,32 

32 I F ( W Z ) 33 , 34 ,27 

33 H=-H 

27 T = ( —ZO + H ) /WZ 
IF(A*T+B >35,34,34 

33 IF (QUADRA <T > >34,99,99 

34 I F ( A*SM+B >36,3 7,37 
37 SM=-(B+SQRTP(E> )/A 

I FLG= 1 
GO TO 99 

36 I F ( QUADRA { SM ) >37,99,99 

29 I F { B ) 22 9 , 99 ,99 

229 IF(-C-EPSI >13,13,97 

12 IF(-E*-*EPSI >13,13,97 

13 I FLG=1 
40 SM = 0 • 0 

GO TO 99 
99 CONTINUE 

50 CONTINUE 

51 CONTINUE 
RETURN 

97 WRITE(6,98) 

98 FORMAT ( 20H CYLINDER SNAFU* ) 

GO TO 51 

END 

$ I BF TC UPDIST LIST 

SUBROUTINE XOYOZO ( DI ST , WX , WY, WZ , XC AP , YCA P, ZCAP, XO , YO , ZO , HZ , I SET , I N 
1 RUP > 

COMMON/ RDE P/LRN 
X00=X0+WX*DI ST 
XTEST=ABS( XOO )-XCAP 
IF(XTEST >910,910,915 
915 DIST=(ABS ( XOO-XO )-XTEST ) /ABS(WX ) 

I NRUP-3 

910 YOO=YO+WY*DiST 

YTEST=ABS (Y00>-YCAP 
I F ( YTEST >920,920,925 
925 DIST=(ABS (YOO-YO) -YTEST ) /ABS( WY) 

I NRUP = 3 

920 ZOO=ZO+WZ*DIST 

ZTEST=ABS( ZOO )-ZCAP 
IF( ZTEST >930,930,935 



93 5 DIST=(ABS ( ZOO-ZO ) -ZTE ST ) /ABS< WZ ) 

I NRU P = 3 
930 CONTINUE 

GO TO ( 940,940,950) , INRUP 
940 X0=X00 
Y0=Y00 
ZO=ZOO 

I F ( ISET.LT.O) ZO=HZ 
950 RETURN 
END 

$ I B FTC ZAP LIST 

SUBROUTINE I NT ACT ( S I GS , S I GT , INRUP) 

common/rdep/lrn 

CALL RANOM(Rlltl) 

IF(R11-(SIGS/SIGT> >810,810,830 
830 I NRUP =4 
810 CONTINUE 
RETURN 
END 

$ I B FTC WHABS LIST 

SUBROUTINE ISOPIK ( N I SO , I I , D , T I , S I GT ) 
COMMON D(8),TI(8) 

COMMON/ ROE P/LRN 
TCS=0.0 

CALL RANDM ( R41 1 1 ) 

DO 4003 NIS=1 ,NISO 
TCS=TCS+D<NIS)*TI(NIS) 

I F 1 R41 -(TCS/S1GT ) >4005,4005,4003 
4003 CONTINUE 
4005 I I =N I S 
RETURN 
END 



Single Level Breit-Wigner Program BWSL 


$ I B JOB 

$ I BFTC SECTBS 

DIMENSION ER( 500),GAMN( 500),GAMG( 500) t GAMF( 500 ) , GAMNO ( 500) 
DIMENSION G ( 100 ) , ER1 ( 100 ) * G AMN1 ( 100 ) , GAMG1 (100) , GAMF1 ( 100 ) , JEND A 
$(100) 

DIMENSION E (2000) ,SIGS( 2000) , SIGC< 2000) ,SIGF( 2000) ,SIGP(2000) , 
$SIGT(2000) , TEMP A (2000), TEMPB( 2000) 

DIMENSION AA ( 10 ) , ABB ( 10 ) ,RR ( 10) ,NOGS ( 10 ) 

DIMENSION CSINT ( 1 0 ) , FS I NT ( 1 0 ) , TOTS I N ( 10 ) 

COMPLEX OEN»EXZ tTX » T2tCEXPl f SUMJT1 f VAL3 
REAL LAM,LAMR 

DATA PI/3. 1415926/, CONS/2. 86E-9/ 

NAMELIST /IN/ NOGS , A A , ABB , RR , I E ND 
READ ( 5 , IN ) 

READ (5,101) LEND, (E ( I ), 1=1, LEND) 

101 FORMAT ( I 6 / ( 1 P6E 1 2 • 5 ) ) 

SIGP(l) = 0. 

SIGS(l) = 0. 

SIGC ( 1 ) = 0. 

SIGFt 1 ) = 0. 

SIGT(l) = 0. 

DO 102 1 = 1,1 END 
A = AA ( I ) 

AB = ABB ( I ) 

R = RR ( I ) 

NGS = NOGS ( I ) 

WRITE (6, 10) I END , I ,AB,A,R 

10 FORMAT! 1H 1 , //////53X , 13, 2X,8H ISOTOP ES/ /////// 30X , 7H ISOTOPE, 10X, 

$9HABUNDANCE,10X,13HATOMIC WE IGHT , 1 OX , 1 2HNUCL E AK KAD . /3 2X , I 3 , 1 OX , 
$1PE14.7,7X,E14.7,8X,E14.7) 

FA = ( (A+1.008986)/A)**2 
DO 100 L = 1 , L END 
NSUMG = 0 
GSUM = 0. 

C SUM = 0. 

F SUM = 0. 

LAM = CONS/SORT ( E ( L ) ) 

BLAM = LAM/PI/2. 

PSIG = FA*AB*PI*1.E24 

JS = 1 

JN = 0 

JE = 0 

JG = 0 

NJ = 0 

J ENDA ( N J ) = 0 

IF (L .GE. 2) GO T012 

JR = 0 

1000 READ (5 ,1) GltJENDf (ERI ( JJ)*GAMG1 ( JJ) » GAMN1 ( JJ) , GAMF1 ( JJ) , J J = 1 , J END 
$) 

1 FORMAT (E12.0,I6/ (4E12.8) ) 

JG = JG + 1 

G ( J G ) = G 1 

JN = 1 + JN 

JENDA (JN) = JEND 

DO 8 J1 = 1 , JEND 

JR = 1 + JR 

I F ( J 1 .EQ. 1 ) JRS = JR 

I F ( J 1 .EQ. JEND) JRE = JR 

ER( JR) = ERHJ1 ) 

GAMN( JR) = GAMN1 ( J 1 ) 



GAMNO(JR) = 0. 

GAMG( JR) = GAMG1 ( J1 ) 

8 GAMF(JR) = GAMFl(Jl) 

WRITE (6, 15 )G1, JEND, ( JJ,ER( JJ ) , G A M N ( J J ) , GAMNO ( J J ) ,GAMG { J J ) , GAMF { J J 
$) , J J- J RS , JRE ) 

15 FORMAT { 1HL, ////////// 46 X,7HG VALUE, 12X , 17HT0TAL MO. OF RES./42X,1P 

$E14.7,15X, I3/////16X,8HRES. MO. , 1 OX , 8HENERG Y 0 , 1 OX, 7HGAMMA N,10X, 
58HGAMMA NO , 10X , 7HG AMMA G , 14X , 7HGAMMA F/ ( 19X , 13,9X,1PE14.7,5X, 
$E14.7,4X,E14.7,5X,E14.7»5X,£14.7> ) 

GO TO 16 

12 JG = JG + 1 
JN = JN + 1 

16 SUM JT1 = (0.,0. ) 

SUMCT1 « 0. 

SUMFT1 = 0. 

JS * JS + JENDA(JN-l) 

JE = JE + JENDA(JM) 

DO 20 J = J S , J E 
L AMR = CONS/SORT (ABS(ER(J) ) ) 

BLAMR = L AMR/ P I / 2 • 

PI = E ( L ) - ER ( J ) 

GAM = GAMN(J) + GAMG(J) + GAMF(J) 

P2 = GAM/2- 

DEN = CMPLXIP1,P2> 

CPF = B L A M R * GAMN< J ) / < Pl**2 + P 2**2) 

T 1 =<GAMN( J >/2./DEN)*BLAMR 
DEBUG LAMR, BLAMR, PI , GAM, DEM, T1 
CT1 = CPF *GAMG ( J ) 

FT 1 = CPF *GAMF(J) 

SUMFT1 = FT 1 +SUMFT 1 
SUMCT1 = CT1 + SUMCT1 
20 SUM J T 1 s T1 + SUMJT1 

DEL = -R/8LAM 
VAL2 = SIN(DEL) 

E X 1 = C M P L X ( 0 * , — D E L ) 

CEXP1 = C EXP ( EX 1 ) 

T 2 = CEXP1*VAL2*BLAM 
ABT2 = C ABS ( T 2 ) 

VAL3 = SUMJT1-T2 
VAL 1 = CABS(VAL3) 

SIGSP = VAL 1 **2 - ABT 2**2 
DEBUG SUMJT1 ,DEL,VAL2,EX1,CEXP1 
DEBUG T2,ABT2,VAL3,VAL1, SIGSP 
94 NSUMG = NSUMG + 1 

GPIECE = G < JG ) *S I GSP 
GSUM = GPIECE + G SUM 
CPIECE = G ( JG ) ^SUMCT 1 
C SUM = CPIECE + C SUM 
FPIEC6 = G(JG)*SUMFT1 

DEBUG I ,L ♦ JG,G( JG) , GP I EC E , G SUM , NGS , NSUMG 
FSUM = FP1ECE + FSUM 
IF (L .GE. 2) GO TO 98 
I F ( MGS .NE. NSUMG) GO TO 1000 
GO TO 99 

98 I F( MGS .NE. NSUMG) GO TO 12 

99 SIG1 = 4.*PSIG*GSUM 
SIG2 = 4.*PSIG*A3T2**2 
SIGP(L) = SIG2+SIGPIL) 

SIGS(L) « S IG1+S IG2+SIGS(L ) 

TEMPA(L) = PSIG*CSUM*BLAM 
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SIGC(L) = TEMP A ( L ) + SIGC(L) 

TEMPB(L) = PSIG*FSUM*BLAM 
SIGF(L) = TEMPB(L) + SIGF(L) 

SIGF(L) = P5IG*FSUM*BLAM + SIGF(L) 

100 SIGT(L) * SIGS(L>+$IGC(L)+$IGF(L> 

CALL ABSIN ( E ,T EMP A , TEMPB, L END ♦ CABS I N, FABS IN) 

CS I NT ( I ) = CABSIN/AB 
FS INT ( I ) = FABS I N/AB 
TOTSIN(I) = CSINT(I) + FSINT(I) 

102 CONTINUE 

WR I TE ( 6 ♦ 1 59 ) 

159 FORMAT C 1 H 1 ) 

WRITE (6,160 ) (I ,CSINT ( I ) ,FSINT( I ) , T0TSIN( I ) , 1 = 1, I END) 

160 FORMAT ( 1H3,27X,64HI SOTOPE CAPTURE FISSION 

1 ABSORPT I0N/2BX , 62HNUMBER INTEGRAL INTEGRA 

2L INTEGRAL/// ( I 3 1 , 1 P E24. 5 , 1 P El 7 . 5 , 1 P El 9 . 5 ) ) 

WRITE(6,80)(LL,E(LL),SIGS(LL),SIGP(LL),SIGC(LL),SIGF(LL) ,SIGT(LL ) , 
$LL = 1 , LEND ) 

80 FORMAT ( 1H1////5X,6HLEVEL, 1 2X , 6H ENERGY , 12X, 7HS I GMA S, 12 X, 7HSIGMA P, 
$ 1 2 X , 7HS I GM A C , 1 2X , ?HS I GMA F , 1 2X , 7HS I GMA T / ( 5 X , I 3 , 1 0 X , 1 P E 1 4 . 7 , 5 X , E 1 
$4.7,5X,E14.7,5X,E14.7,5X,E14.7,5X,E14.7i ) 

WRITE (6,111 ) (E( I ) , S I GC ( I ),SIGS( I),SIGF( I ), 1 = 1, LEND) 

111 FORMAT( 1H$,1P4E12.5) 

STOP 

END 

$ I B FTC A B S INI LIST 

SUBROUTINE ABSIN ( E , TEMP A ,TEMPB ,LEND, CABS I N, FABS IN) 

DIMENSION E ( 2000 ), TEMP A ( 2000 ), TEMPB ( 2000 ) 

CABSIN = 0,0 

FABSIN = 0.0 

DO 100 1=2, LEND 

XNAME = ALOG ( E (I > / E ( I — 1 ) ) 

CABSIN = CABSIN +((TEMPA(I) + TEMP A ( I -1 ) ) / 2 . )*XNAME 
IOC FABSIN = FABSIN +((TEMPB(I) + T E MPB ( I -1 ) ) / 2 . ) *XNAM E 
RETURN 
END 

SDATA 
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Doppler Broadening Program DBCS 


$ I B JOB 

SIBFTC DBCS LIST 

COMMON NLEV, I , J,Z,CON,DBI ,E,SIGC,SIGS,SIGF, ITYPES,NEDU,NEDL 
COMMON EPSUtKSAV 

DIMENSION E ( 2 200 ) , S I GC ( 2200 ) , S I GS ( 2200 ) , S I GF ( 2200 ) , S I GT ( 2200 > , 

$DB IC ( 2200 ) , DB I S ( 2200 ) ,DBIF(2200) 

EXTERNAL DBO 

1000 RE AD ( 5 » 1 01 0 ) A , T , T END , NL EV , I T YP ES t NPUNCH , NEDU , NEDL , I NT ER 
1010 FORMAT! 1P3E 12.5, 6 16) 

READ! 5,1020) < E ( I ) , S I GS ! I ) , S I GT ! I ) ,1=1, NLEV) 

1020 FORMAT t 1P3E12.5 ) 

DO 10 JK = 1 ,NLEV 
10 SIGC! JK)=SIGT(JK)-SIGS< JK ) 

KSA V = 1 
JINTER = NEDU 
CALL TIMEl(TS) 

Z = A/ ( 8.6167E-5*T ) 

CON= SORT ( A/ ! 1.0828E-3*T) ) 

DO 2000 I =NEDU , NE DL 
1=1 

I F ! I .NE. JINTER) GO TO 1050 
JINTER = JINTER + INTER 
CALL T I ME 1 ! T P ) 

TELAP = ! T P-T S )/36 00. 

I F ( TEL AP .LT. TEND) GO TO 1050 
NEDL = J INTER-INTER-1 
GO TO 2010 

1050 EPSL= ! < SORT (E(I))-SORT (88.0/Z))**2) 

EPSU=( (SORT (Eli) J+SORT (88.0/Z) ) **2 ) 

IF(EPSL— E(NLEV) ) 1540, 1550, 1550 
1540 E PSL = E ( NL EV ) 

WRITE (6,3050) E(NLEV),I 

3050 FORMAT! 1 HL , 1 OX , 49HTHE LOWER INTEGRATION LIMIT HAS BEEN SET EQUAL T 
101 PE 12. 5 ,12H FOR SET NO. 16) 

1550 CONTINUE 

IF ( EPSU-E ( 1 ) ) 1570, 15 70, 1560 
1560 E PSU = E ( 1 ) 

WRITE(6,30 75) E(1 ),I 

3075 FORMAT ( 1HL,10X,49HTHE UPPER INTEGRATION LIMIT HAS BEEN SET EQUAL T 
1 01 PE 1 2 • 5 , 1 2H FOR SET NO. 16) 

1570 CONTINUE 

NTYPES-ITYPES+2 
DO 1090 J = 1 ,NT YPES 
J = J 

DB I =E XPS 1 (EPSL,EPSU,DBO,KK ) 

GO TO ( 1060,1070,1080) ,J 
1060 DB IC ( I ) =DB I 
GO TO 1090 
1070 DB I S ( I ) =DB I 
GO TO 1090 
1080 DB I F ( I ) =DB I 
1090 CONTINUE 

S IGT ( I ) ~DB I C ( I ) +DB I S ( I ) +DB I F ( I ) 

2000 CONTINUE 
2010 WR I TE ( 6 , 3000 ) NEDL 

3000 FORMAT ( 1 HI //////40X , 16HN0. OF ENERGY S =,I5 ) 

WRITE(6,3005) A , T 

3005 FORMAT ( 1HL »5X ,43HTHESE ARE CROSS-SECTIONS FOR ISOTOPE NUMBER1 PE 1 2 • 
15/1H ,2X,20HD0PPLER BROADENED TO 1 P E 1 2 . 5 , 1 6H DEGREES KEL V I N . // / / 1H 
2 ♦ 3X , 6HENERGY ,9X, 7HS IGMA C,8X,7HSIGMA S,8X,7HSIGMA F,8X,7HSIGMA T) 



WRITE (6,3010 ) ( E ( I ) ,DBIC( I ) ,DBIS(I> ,DB1F( I ) ,SIGT ( I ) * I=NEDU,N£DL ) 
3010 FORMAT ( 1P5E15.5 ) 

I F ( N PUNCH )30 60 *3030*3020 

3020 WRITE ( 6,302 5 ) ( E ( I ) , DB I S ( I ) , S I GT ( I ) , I=NEDU, NEDL ) 

3025 FORMAT ( 1 H$ , 1 P3E 1 2 • 5 ) 

GO TO 3030 

3060 CALL BCDUMP (E (NEDU) ,E (NEDL ) ) 

CALL BCDUMP ( DB IS (MEDU)»DBIS( NEDL) ) 

CALL BCDUMP (SIGT( NEDU) , SIGT(NEDL) ) 

3030 CONTINUE 

GO TO 1000 
END 

$ I BFTC EXPS1. LIST 

FUNCTION EXPS1 (XMIN»XMAX,FUNC1»KER) 

DIMENSION V (500 ) ,H(500 ) , A( 500 ) » B ( 500) t C( 500) ,P( 500 ) , E ( 500 ) f NE ( 500 ) 
EQUIVALENCE ( E,NE ) , (TEST * NTE ST ) 

T=3 .0E-5 
V( 1 ) = XM I N 

Hi 1 ) =0 • 5* ( X M AX“XM IN) 

C< 1 ) =FUNC1 ( XMAX ) 

B ( 1 ) = FUNC 1 (XMIN+H( 1 ) ) 

A ( 1 ) = FUNC 1 ( XMIN ) 

P(1)=H(1)*(A( 1 )+4.0*B( 1 )+C( 1 ) ) 

E ( 1 ) =P ( 1 ) 

A N S = P ( 1 ) 

N=1 

FRAC=2.0*T 

1 FRAC=0.5*FRAC 

2 TE ST=ABS ( FRAC*ANS ) 

K = N 

3 DO 7 I = 1 » K 

4 IF(NTEST-IABS (NE ( I ) ) )5,5,7 

5 N=N+1 

V( N)=V( I ) +H ( I ) 

H(N )=0.5*H( I ) 

A ( N ) = B ( I ) 

B(N)=FUNC1 <V(N)+H(N) ) 

C(N)=C( I ) 

P ( N ) =H ( N ) * ( A ( N ) +4 #0*B ( N )+C ( N) ) 

Q=P ( I ) 

H(I)=H(N) 

B ( I ) =FUNC1 ( V ( I ) +H ( I ) ) 

C ( I )=A (N) 

P(I)=H(I)#(A(I )+4.0*B ( I ) +C ( I ) ) 

Q = P ( I ) +P ( N ) — Q 
ANS= ANS+Q 
E ( I ) =0 
E ( N ) =Q 

6 IF(N— 500)7 ,13,13 

7 CONTINUE 

8 I F ( N— K ) 9 ,9 , 2 

9 0=0.0 

10 DO 11 1=1, N 

1 1 Q=Q+E ( I ) 

12 IF<ABS(0)-T*ABS(ANS) ) 14,14,1 

13 KER=KER+1 

14 ANS=0.0 

15 DO 16 1=1, N 

16 AN S = ANS + P ( I ) 

EXPS1=( ANS+0/30.0)/3.0 



17 RETURN 
END 

$ I BFT C DBO * LIST 

FUNCTION DBO (XX) 

COMMON NLEV, I , J,Z ,CON,DBI ,E,SIGC,SIGS,SIGF, I T Y P ES r NEDU , NEDL 
COMMON XMAXtKSAV 

DIMENSION E ( 2 200 ) ,SIGC(2200) ,SIGS(2200) ,SIGF(2200> 

X = XX 

DO 2300 K =KS AV , NL EV 
IF( E (K)-X >2050,2300, 2300 
2050 GO TO ( 2060,2070, 2080 ), J 
2060 SIGM=SIGC(K) 

S I GN=S I GC ( K— 1 ) 

GO TO 2100 
2070 S I GM=S I GS ( K ) 

S I GN=S I GS ( K- 1 ) 

GO TO 2100 
2080 S I GM=S I GF ( K ) 

S I GN=S I GF ( K- 1 ) 

2100 CONTINUE 

EM= ( S IGM-S IGN ) / ( E ( K ) -E ( K—l ) ) 

B=S IGM-EM*E ( K ) 

SIGX=EM*X+B 

ARG1 Z* ( (SORT ( E < I ) )-SQRT (X ) )**2) 

ARG2=-Z*( (SORT (E( I) ) + SQRT (X) )**2> 

PI = EXP ( ARG1 ) 

P2 = E X P ( ARG2 ) 

D BQ= ( ( CON/E ( I) )*SGRT(X)*SIGX)*(EXP(ARG1)-EXP(ARG2) ) 

IF ( X • E 0 * XMAX) K S AV = K — 1 
GO TO 2350 
2300 CONTINUE 
2350 CONTINUE 
RETURN 
END 

$DAT A 
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TABLE I. - DOPPLER COEFFICIENTS FOR GOLD 197 


Square root 


Doppler coefficients, a 


Effective 

of surface to 





resonance 

mass ratio, 

Monte Carlo 

Experimental 

ZUT 

TRIX 

integral, 

Vs/M 





T eff 

1.30 

a 4. 06X10 -5 

b 3. 15X10" 5 

4. 31X10' 5 

b 3.49xl0' 5 

130 



4 95 


3 95 


2.85 

3.99X10' 5 

b 3. 21 

3. 69X10' 5 

b 3. 41 

270 



3.45 


2. 90 


6.02 

6. 34X10’ 5 

5.35X10" 5 

5.00 

550 

6. 99 

6.93 


6.08 

5.35 

610 

9.89 

7.06 


6.71 

5.60 

800 

11.4 

6.93 

b 4. 62X10 -5 

6.65 

5. 62 

868 

14.0 

5.69 

4.70 


4.80 

990 

25.4 

C 4. 27 


3.26X10" 5 

2.75 

1286 

42.5 

C l. 98 

3.75X10' 5 

1.61 

1.30 

1465 


a This Monte Carlo value is corrected to allow for an upper cutoff of 500 eV 
(J) used. The value of AI^ to this energy was 3. 13, which represents 
only 77 percent of the total. 

b These values are from ref. 11. All other experimental and TRIX values 
are from ref. 12. 

c The effect of different cross-section sets was observed by rerunning each 
case with the ZUT agreement (see table n). 


TABLE H. - MONTE CARLO AND ZUT RESULTS FOR IDENTICAL 


CROSS-SECTION SETS 


Square root 

Monte Carlo 

ZUT 

Difference, 

percent 

of surface to 
mass ratio, 

V S/M 

Change in 
effective 

Effective 

resonance 

Change in 
effective 

Effective 

resonance 

resonance 

integral, 

resonance 

integral, 



integral, 

I eff 

integral, 

Z eff 



AI eff 


AI eff 



25.4 

10.05 


10.3 


2.4 

1163.43 

1183.7 




1.7 

42.5 

5. 64 

5.72 



1.35 

1321.72 

1331. 67 




.75 
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Doppler coefficient, a D (°K) 


i 
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Square root of surface to mass ratio, Vs/M, (cm)(g _1/2 ) 


Figure 1. - Calculated and experimental values of Doppler coefficients as function of square root of 
surface to mass ratio. 
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shall provide for the widest practicable and appropriate dissemination 
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